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Abstract 

We consider the product of M quadratic random matrices with complex elements and no further 
symmetry, where all matrix elements of each factor have a Gaussian distribution. This generalises 
the classical Wishart-Laguerrc Gaussian Unitary Ensemble with M = 1. In this paper we first 
compute the joint probability distribution for the singular values of the product matrix when the 
matrix size and the number M are fixed but arbitrary. This leads to a determinantal point process 
which can be realised in two different ways. First, it can be written as a one-matrix singular value 
model with a non-standard Jacobian, or second, for M > 2, as a two-matrix singular value model 
with a set of auxiliary singular values and a weight proportional to the Meijer G-function. For 
both formulations we determine all singular value correlation functions in terms of the kernels of 
biorthogonal polynomials which we explicitly construct. They are given in terms of hypergeometric 
and Meijer G- functions, generalising the Laguerre polynomials for M ~ 1. Our investigation was 
motivated from applications in telecommunication of multi-layered scattering MIMO channels. We 
present the crgodic mutual information for finite- A^ for such a channel model with M — 1 layers of 
scatterers as an example. 
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1 Introduction 



Random matrix theory (RMT) remains a very active field of research, after many decades of work. 
While originally being conceived in the area of mathematical statistics and nuclear physics, today's 
applications of RMT extend beyond the mathematical and physical sciences even in a broad sense, 
and we refer to [1] for a recent overview. 

One of the topics in RMT that has caught recent attention is that of products of random matrices. 
Having one of its original motivations in statistical physics in the description of chaotic and disordered 
systems 0, among more recent applications are combinatorics [3] and telecommunications [1], which 
has also been one of our motivations. In particular, we consider MIMO (multiple-input and multiple- 
output) communication networks, where multi-antenna transceivers are utilised to improve the system 
capacity in a rich scattering environment. 

Products of matrices loose much of the symmetry of the individual matrices and are generically 
complex. For simplicity we will consider the individual matrices to be complex, too, with independent 
Gaussian distributions. The spectral properties of matrix ensembles carry important information. It 
is encoded in the eigenvalue decomposition as well as in the singular value decomposition. We will 
focus on the latter here. 

A striking property of RMT is its universality, that is the independence of the underlying dis- 
tribution of the individual matrix elements. It is usually manifest in the limit of large matrix size. 
However, if we study the local, microscopic behaviour of the spectrum on the scale of the mean level 
spacing between singular values, it is often vital to have a detailed knowledge of the joint distribution 
of singular values (or eigenvalues) at hand for finite matrix size. In particular to derive a determi- 
nantal or Pfafhan structure of the correlation functions of the random matrix ensemble has proven 
very useful for universality studies. Some of the most powerful proofs of universality start from the 
knowledge of orthogonal polynomials of these determinantal or Pfaffian point processes, in order to 
perform the asymptotic analysis. We refer to [5] for a standard reference. 

The aim of this work is to provide such a starting point, by deriving the joint probability density 
function (jpdf) of the singular values of the product matrix at finite matrix size N, for a finite product 
of M matrices. We consider the simplest case of M quadratic N x N matrices of Wishart type, 
independently and identically distributed by Gaussians with unit variance. In telecommunications 
this is also the setting often encountered, with both N and M finite. The singular value distribution 
of products of complex Wishart matrices is then the setup for the calculation of several information- 
theoretic quantities. 

In previous works the spectral density of singular values as well as the moments of such product 
matrices were derived in the macroscopic large- limit. They use probabilistic methods such as free 
random variables O [71 [8], field theoretic methods such as planar diagrams |9], and inverse Mellin 
transforms [3]. The limit of infinitely many matrices in a product were studied in other works, either 
for finite-size \W\ [TT| \T7[ [T3] or for infinitely large matrices \14:\ [T5] , where the problem was mapped 
to a differential equation. 

Very recently the jpdf and its correlation functions of the complex eigenvalues of the matrix 
ensemble we are considering has been derived for finite N and M \16 \ \n \ \T8\. Also here the macroscopic 
large- A^ density in the complex plane was known previously, see |19| [20l [9] for a collection of works. 
However the corresponding question about the singular value correlation functions for finite N and 
M was still open and is addressed in our work. 

This article is organised as follows. In section [2] we determine the jpdf of singular values, using 
two different ways. Section [3] is devoted to the computation of the correlation functions, by first 
determining the biorthogonal polynomials associated to this problem in subsection 13.11 and then the 
kernel(s) leading to all fc-point correlators in subsection 13.21 The spectral density itself is discussed 
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in more detail in subsection 13.31 where we compute all its moments and identify the correct rescaling 
for the macroscopic large- limit of our results. In particular we compare our results to the known 
large- asymptotic spectral density e.g. from [9]. Section |4] illustrates how our results can help in 
applications in telecommunications, by computing the ergodic mutual information and comparing it 
to numerical simulations. After presenting our concluding remarks and open questions in section [5] we 
collect some technical tools in two appendices. 



2 Joint probability distribution of singular values 

We are interested in the singular values of the product Pm of M independent matrices with complex 
matrix elements of size N x N, Xj G G1(A^, C) for all j = 1, ... , M: 

Pm = XmXm^i ■ ■ ■ Xi . (2.1) 

Furthermore the matrices Xj are distributed by identical, independent Gaussians, 

V{Xj) = exp [-TVXjATjj . (2.2) 
The partition function Z^^'^ is then defined as 



M 

(M) 



C I f[d[X,]P{X,) , (2.3) 



where d[Xj] = n^/3=i d{Xj)a/3d{X*)a/3 denotes the flat measure over all independent matrix elements. 
In this section we compute the jpdf of the singular values of Pm- For M = 1 this is the well known 
Wishart-Laguerre (also called chiral Gaussian) Unitary Ensemble, and for most of the following we 
will thus restrict ourselves to M > 1. We will not keep track of the normalisation constant C here, 
and will only specify it later, once we have changed to singular values. 

In the first step we perform the following successive change of variables from Xj to Yj : 

Yi = Xi, and Yj = XjYj^i for j = 2, . . . , M , (2.4) 

e.g. Y2 = X2X1, I3 = X^X2Xi etc. In the new variables Ym = Pm is the product matrix we are aiming 
at. While the first one, Yi = Xi, is a trivial relabelling, each subsequent change of variables carries a 
non-trivial Jacobian given by 1/ det[y/_iyj_i]^ for j = 2, . . . , M. This can be seen as follows. Due to 

d{Yj)ai3 = Z^^i d{Xj)a'y{Yj^i)jp every column vector of the matrix Yj acquires a factor 1/ det[Yj_i] 
from the change of variables, and likewise its complex conjugate. Taking into account all A^ column 
vectors and their complex conjugates we find the given Jacobian for each j = 2, . . . , M, and we arrive 
at 

„ M M 

4*'^ = CJ exp [-TtYIy,] J] _^__^exp [-Tr i;ty^.(i^^^y^._,)-i] . (2.5) 



i=l 



In writing this we assume that the matrices Xj and thus their products Yj are invertibl^. In the second 
step we decompose each matrix Yj = VjAjUj, j = 1, . . . , M in its angles and singular values, where 
Aj = diag(A[''\ . . . , A^^) contains the positive singular values G M+, a = 1, . . . , A^, and Uj G U(A^) 



^Our restriction from general complex N x N matrices to Gl(Af, C) removes only a set of measure zero. 
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and Vj £ U(A'^)/U(1)^ are unitary. The Jacobian resulting from the singular value decomposition of 
each matrix Yj is well known and is given in terms of the Vandermonde determinant, 



Af>a>6>l 



b > 



det 

l<a,h<N L 



(2.6) 



Since we encounter the matrices Yj in the combination Y^Yj, only, the unitary matrices Vj completely 
drop out which leads to 



Z 



(M) 
N 



M 



N 



C d[VM[U^] n \ exp [-TY A?] A^(A? 

1 



2\2 



i=l 



a=l 



M 

» A/ r N 
C Yl\ dmd[U^] n d^'a^^^^ \ [-TV A?] A^(A? 

1=1 { a=l 



2\2 



2\2 



M 

i=2 



det[4_J^ 



-Tr C/j A2C/,.A7_\ 



(2.7) 



In the second step we employed the invariance of the Haar measure d[Uj] under Uj — )■ UjUj-i, which 
leads to the decoupling of the integrations over d[Ui] and all the d[Vj] from the rest of the integrals. 
The remaining unitary integrations in the last line of eq^ (|2.7p can be performed using the so-called 
Harish-Chandra-Itzykson-Zuber (HCIZ) integral |2H | 



d[Uj]exp -Tr f f/jA^f/jA-^i 



A7v(A2)A7v(AT^J l<a,b<N 



exp 



{21 



The Vandermonde determinant of inverse powers is proportional to the ordinary one with positive 
powers due to the following identity: 



Ajv(A: 



det 

l<afi<N 



0')\2fe-2 



det[A2]^-i 



(2.9) 



This leads to many cancellations in eq. (12. TP . in particular of almost all Vandermonde determinants: 



Z 



(M) 
N 



•^0 a=l [ .7=1 Aa 

M 

X TT det 



(i) 





AT 




1 exp 




A;v(A?)A;v(Ai^ 




. 6=1 





1=2 



l<c,d<Ar 



(A^' 



(2.10) 



We expand the determinant comprising A^^^ and Ai^'' in N\ terms. Each of these terms yields the same 



(aJ-^^. 

an 

contribution since the permutation involved in the definition of the determinant can be absorbed in the 

(2) (3) 

determinant comprising A^ and Ac due to the antisymmetry of determinants, and a relabelling of the 

(2) (3) 

integration variables. Next we expand the determinant comprising A^ and Ac whose permutations 



The Haar measure is normalised such that there is no further proportionahty constant on the right hand side. 
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can be absorbed in the determinant comprising A^^'* and and so on. This interplay of expansion 
and absorption of the permutations of the determinants can be continued untih all determinants 
stemming from the HCIZ integral are replaced by their diagonal part. Note that we do not require 
any symmetrisation in the variables aI*^^ here. Hence our argument applies to correlation functions 
of the aI^^^ in the next section [Sj too, where we integrate the joint probability distribution function 
(jpdf) only over a subset of these singular values. 

Almost all remaining multiple integrals can be simplified as follows: 



N 

n 

a=l 



Ai^) 



exp 




exp 



U)\2 



X det 

l<c,d<7V 



det 

l<c4<N 



det 

l<c.d<N 



)\2d-2 



dX. 



(1) 



(Ai^^)3-2'i 



1 r<M, 



exp 



2M-1^0,M I 0,...,0,d-l 




dX 



A 



(i) 



exp 



(A 



exp 



(i)^2 



{X'a 



(A^)^2 



ix^'-"^) 



(A 



exp 



(A^ 



(A/-l).2 



(A 



(2.11) 



Notice that we have left out the integration over the variables aI^^^ . In the second line of eq. (12. lip we 
have pulled all the integrations into the corresponding rows of the determinants, and in the third line 
we have used the integral identity (jA.lOh in the squared singular values, that is derived in appendix Rl 
The special function appearing here is the so-called Meijer G-function, see eq. ()A.ip for its definition 
|23j . The number of zeros in the bottom line of the Meijer G-function is M — 1. Our first main result 
is thus the following singular value representation of the partition function, after changing to squared 
singular values Sa = (Ai*^-*)^ with dsa = 2X\^^'' dX^a '^\ a = 1, . . . , A^, in eq. ()2.10p : 



7(m 

■'N 



^jpdf(s) 



CF^ T\dsa^N{s) det 

n l<c,d<N 
a=l - ' - 



M,0 

0,M I 0,...,0,d-l 



N 



ci^'^AMis) det 



^M, 

^0,M \ 0,...,0,ci-l 



n d^aVjpdiis) , (2.12) 
(2.13) 



a=l 



where "Pjpdf is the jpdf. We will show later that it corresponds to a determinantal point process. The 
constant in front of eq. (|2.13p . 



N 



\M+1 



(2.14) 



a=l 



has been chosen such that the partition function is normalised to unity. This can be seen as follows. 
Applying the Andreief integral identity. 



N 

y n J' 



det [(/.e(sd)] det [V'c(sd)] = A^! det 

d<N l<c4<N l<c4<N 



ds 4>c{s)i)d{s) 



(2.15) 



which applies to any two sets of functions (j)c and ijjc such that all integrals exist, we obtain for the 
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partition function 



rim 

-'N 



'N ^ 



C7^"^iV! det 

l<c4<N 



Ubb Ltq, A/ I 0,...,0,d-l 



C^'^iV! det \r{cf^-^r{c + d-i)] 

l<c,d<N '- 

Af N 



(2.16) 



c=l 



6=1 



Here we have used another integral identity from the appendix, eq. ()A.7p . and pulled out factors of 
Gamma- functions from the rows of the determinant. The remaining determinant is nothing but the 
normalisation of the Wishart-Laguerre ensemble (at M = 1) which is well-known |24j . 

We refer to the result (j2.12p as a "one-matrix" singular value model because it is of the form 
that would result from the singular value decomposition of a single random matrix, however with a 
non-standard Jacobian ^ Ajv(s)^. As an easy check we can see that due to 



1,0 

0,1 I d-l 



„d-l 



exp[ 



(2.17) 



the expression in eq. (12.12P reduces to the standard Wishart-Laguerre ensemble when setting M = 1, 
and taking the exponentials out of the second determinant. In principle we could now try to compute 
all singular value fc-point correlation functions defined as 



(N-k) 



■ a=k+l 



(2.18) 



For example for k 
following 



1 this gives the spectral density which is normalised to in our convention 



N 



(2.19) 



whereas for A; = A we have the jpdf itself, R)^ {si, . . . , sn) = N\Vjpd{{s)- However, due to the matrix 
inside the second determinant in eq. (j2.12p being labelled by indices of the Meijer G-function, the 
computation of the R^^ is a highly nontrivial task. We postpone this computation to section [3] using 
a second "two-matrix" formulation, that is introduced in the next subsection. 



2.1 An alternative jpdf with auxiliary variables 

We introduce a formalism that is more convenient to handle when computing correlation functions. 
Let us step back by considering eq. (I2.10p . taking for simplicity M = 2: 



^(M=2) 
-'N 



c 



N N ,,(1) 



a=l 



X det 

l<c4<N 



exp 



(A 



b=l H 



-(A, 



{1)n2 
b ) 



A;v(A?)A;v(Al) 



(Al^^)^ 



(2.20) 



This is precisely of the form of a "two-matrix" singular value model (2mm) that results from the 
singular value decomposition of a two-matrix model, see e.g. in [25]. The advantage is that now 
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we have the standard form of the Jacobian given by one Vandermonde per set of variables and an 
additional determinant of a function that couples the two sets of variables. Such a setting can be 
tackled using the known biorthogonal polynomial technique, as reviewed in |26] . In order to apply 
this technique we have to bring eq. (j2.10p into such a form, but for arbitrary values of M > 2. This can 
be readily achieved by taking the same step as from eq. (j2.10p to eq. (j2.1ip . but this time excluding 
both sets of integrations over the variables and aI^^. 

The symmetry argument goes along the same lines replacing determinants by their diagonal parts, 
see discussion after eq. (I2.10p . but this time we keep the determinant containing the exponential with 
A^^^ and A^^^\ The integrations are: 




(xi^^? 



det 

l<c,d<N 



exp 



exp 



(A 



(A 



exp 



det 

l<c4<N 





(aW)A1 


^0,...,0 





( ixf^A] 






(Ai^))^V 









(2.21) 



where we have used again the identity (jA.lOp from appendix [Aj see also [16] for a related recurrence 
relation. We therefore arrive at our second main result, the following 2mm representation for the jpdf, 
after changing again to squared singular values Sa = (aI*^'')^, ta = (aI^'')^, a = 1, . . . , N: 



■'N 



c 



(M) 
N 



N 



\{ds, 



dta 



" a=l 

N 

l[dSadtaVf^Tis,t) 



e"*" AN(s)AN(t) det 

l<c,d<N 



M-1,0 

0, A/-1 \ 0,...,0 



a=l 



c 



(A/) N 



Vl^Tis,t) ^ ^Ht-'e-'- Ar,is)A^it) det 

■"^ iV! -^-^ l<c,d< 



a=l 



d<N 



A/-1,0 

0, A/-1 I 0,...,0 



(2.22) 



It is normalised to unity as we will check below. 

The crucial advantage in comparison to "Pjpdf is the matrix inside the determinant which has indices 
that label the integration variables, and not the indices of the Meijer G-function as in eq. (|2.12p . This 
2mm describes the correlations among the singular values aI^^ of a single matrix Xi and the singular 
values aI^^ of the entire product matrix Pm, to be computed in the next section [3l 

As a check for M = 2 we get back to eq. (j2.20p . using 



1,0 _ 
0,10 



(A 



(2)^2' 



exp 



(A 



(2)n2' 



(2.23) 



Confirming the normalisation in eq. (12.22P is at the same time a check that this representation can 
be mapped back to eq. (j2.12p in a different way. Applying once again the Andreief formula (|2.15p to 
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eq. (I2.22p . but this time only to the t-integration, we obtain the following: 

TV 



/ TT ^e-*'' A^(i) det 

Jo ta ^ \<c,d<N 



G. 



M~1,0 

0, M-1 \ 0,...,0 



Nl det 

l<c4<N 



Nl det 

l<c4<N 



ULL L Ltq, I 0,...,0 



'-'0,M \ 0,...,0,d~l 



(2.24) 



upon using the identity ()A.9P with m = M — 1 from appendix [Aj This brings us back to the "one- 
matrix" model representation eq. (I2.12P in terms of a single set of singular values, with the proper 
normalisation. 



3 Singular value correlation functions 

We are now prepared to compute arbitrary fc-point correlation functions of the singular values. Rather 
than using the definition (j2.18p we will consider the more general correlation functions of the jpdf in 
the 2mm representation (I2.22P : 



^k^H^i^ ■ ■ ■ T Sk;ti, . . . ,ti) = 



{N-ky.iN-iy.j, 



, N N 

n n dhviTi^^t). (3.1) 

a=k+l b=l+l 



The fc-point functions of the singular values of the product matrix Pm in eq. (|2.18p can be obtained 



by integrating out all auxiliary variables, or by setting I = 0: ii[.^Q^(si, . . . , s^; — ) = ^k^'i^i^ • • • > ■Sfc)- 
Let us introduce the following set of biorthogonal polynomials (bOP) in monic normalisation, 
Pi{x) = + . . . and qj{x) = + . . 



dsdtw'^^''){s,t)pf'\s)qf'\t)=5,,hf'\ 



AM) 



with squared norms h^j^^ and the weight function defined as 



w 



0, A/-1 I 0,...,0 



(3.2) 



(3.3) 



for M > 1. These polynomials are guaranteed to exist following the general theory of bOP that was 
very recently further developed [27], see also |26] for a recent review. We will explicitly construct the 
bOP for general M > 1. The general (/c,Z)-point correlation functions eq. (j3.ip are given in terms of 
four kernels that are constructed from the kernel of bOP [281 \29\ I30j : 



^k]P(^^^ • • • jSfc;*!, • • • ,ti) 



det 



Hoi{Sa,Si,) HoQ{Sa,tj 



l<a,b<k 



l<a<k 
l<j<l 



Hii{ti,S}j) Hio{ti,tj) 



l<b<k 
l<i<l 



(3.4) 



The kernel of bOP is defined as 



j=0 



(3.5) 
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All four kernels Hab are based on this relation, 

Hoi{Sa,Sb) 



Hoo{Sa,tj) 



/"OO /"OO / 



ti 



0,M-1 I 0,...,0 



M-1,0 

0, Af-1 I 0,...,0 

OO 



poo 



0,...,0 



(3.6) 
(3.7) 

(3.8) 
(3.9) 



Note that eq. ()3.4p implies that both the one- and two-matrix model jpdf represent determinantal 
point processes, i.e. for the former eq. (j2.13p becomes 



^jpdf(s) = ]V!,<det<^[^0i(sa,56)] . 



3.1 The biorthogonal polynomials 

In order to compute the bOP let us first determine the bimoment matrix 



(3.10) 



ds / dtw^^'\s,t)sH^ 

JO 

JO * 

OO 

dtt^+'e-'{i\)^^-^ 



(i + j)!(i!)^^-i , 



0, A/-1 I 0,...,0 



(3.11) 



which follows again from an identity for Meijer G-functions, see eq. (jA.7P appendix El The bOP as 
well as their norms are determined by this bimoment matrix (see e.g. |27j): 



D 



(A/) 



det 



-^00 

loi 



ho 
hi 



InO 
Inl 



D 



(A/) 



det 



hn—1 ■ ■ ■ Inn—1 

1 S ... 

-^00 ho • • • In-lO 1 

IqI III • • • In-U t 



(3.12) 



where 



DiM) ^ 



'On -'In 



n-1 



'n— In 



(3.13) 



0<J,j<n— 1 -'■-^ 0<ij<n— 1 

r,m/r,{M) 



i=0 



(3.14) 
(3.15) 



9 



In order to have more explicit expressions it is instructive to compare these equations with the standard 
Laguerre polynomials L„,(x). We need them in monic normalisation denoted by L„(x): 



LJx) 



k=0 



{n-ky. \kl 



^ X 



with squared norms 
and a symmetric bimoment matrix 

1 

dtf+^e-' = {i+jy. , K{x) 



dxe ''Ln{x)Lrn{x) = 5nm{n\f = Snm , 



(3.16) 



(3.17) 



^(M=l) 



det 



-^00 



A/=l 



Ion\ 



M=l 



X 



(3.18) 



The former equals eqs. (13. IIP with M = 1. As can be seen the monic Laguerre polynomials also 
have a determinant representation from the Gram-Schmidt procedure, which is exactly the one in eq. 
(13.12|) at M = 1 (or eq. (I3.13P as they become equal then). 

From the comparison of eqs. (|3.12p and (j3.13p to eq. (|3.18p we can read off the following. For 
the qn{t) we take out the common factors from the first n columns of the determinant, 

i = 0, l,...,n — 1, with the remaining determinant being identical to that of the monic Laguerre 
polynomials. The determinant of the bimoment matrix (I3.14P is already written to be proportional to 
the corresponding one of the Laguerre ensemble. We thus have 



^^^^ Lnit) 

p^n-l(,,)A/-l 



Ln{t) , 



(3.19) 



for all values of M. Likewise we can read off the squared norms by comparing them to the Laguerre 
case eq. (I3.17p : 

hW = ffi#)!^/,(M=l) = („,)A^+1 . 



This equation is formally redundant for M = 1. 

For the polynomials Pn{s) the case is slightly more complicated. Also for these polynomials we 
can take out common factors from all n + 1 columns, however this will modify the arguments in the 
last row of the determinant in the numerator: — )• sY(z!)*^~^. Expanding with respect to the last 
row we get a polynomial with the same coefficients as the monic Laguerre polynomials, but now with 
instead of the monomial s* alone, resulting into 



\\\M-1 



(3.20) 



hA/-l 



(-ir 



E (n-k)\ 
(-l)"(n!)^^ iFA/(-n;l, 



fc=0 

k\ ) 



{n-ky. 



(fc!) 



A/-1 



M+1 



hs). 



(3.21) 



This function can be interpreted as a generalisation of the monic Laguerre polynomials reobtained 
when setting M = 1, see eq. (|3.16p . Moreover we could express it in terms of the generalised 
hyper geometric function iF^j which has M arguments equal to 1 in the second set of its indices. 

The kernel (j3.5p is now completely determined. We proceed by computing the various kernels (j3.6p . 
(13. Sp and (13. 9p by integrating the kernel of bOP. 
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3.2 The kernels and all correlation functions 

We start with the kernel Hqi that is relevant for the correlation functions of all singular values 
s„ = (Ai^^)2. We have 



.-n h\ Jo \ 



j=0 "-j 

7V-1 



Sb 



AM) 

where we introduce the following integral transform 



(3.22) 



dt'Lj{t')t'-\-''G^']^lL\ 



o,...,o 



^ (-ly 

^ (j -i)! V^! 

i=0 ^ ^ 

j=0 



0,...,0 



A/, 



1 I „•! ; '^0,M I 0,...,0,i 



S6 



(3.23) 



upon using the identity eq. ()A.9P for d — l = i and m = M — 1. A more compact expression of can 
be found by combining the Rodrigues formula 



LAt') = e' (-^) 



d 



(3.24) 



and the identity (|A.13h . We substitute t' — ^ Sb/t' in eq. ()3.23p and express the derivative in eq. ()3.24p 
as a derivative in Sh- The integration over t' yields 



0,...,0 



Sb 



M, 1 



Af+1 I 0,...,0 



(3.25) 



The second equality can be obtained by applying the definition (jA.ip of Meijer's G-function. We can 
thus write down the full answer for all /c-point correlation functions of the singular values: 



Rf^\si,. . . ,Sk) = Rl%\si,...,Sk;-) = det [i/oi(sa,Sb) 

' l<.a.b<k 



det 

l<a,6<fc 



^ — iFui-j; 1, ■ ■ ■ , 1; Sa)Gf' 



M, 1 / -i 

M+1 1 0,...,0 



Sb 



(3.26) 



The simplest example is the density or 1-point correlation function of singular values which is given 
by 



N-l 



r[^\s) = Hoi{s,s) = Y,\ii'Mi-j;h---A;s)Gl 



A/+1 \ 0,...,0 



(3.27) 



This example will be further discussed in the next subsection 
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The next kernel Hqq is readily given from its definition eq. (13. 7p together with eqs. (I3.2ip . (I3.19P 
and (j3.20p . We therefore turn to Hiq from eq. ()3.9p 



N-l ^ / /• 



as Pi [S)Uq ;^j_^ I Q 



1=0 

Af-1 



Liitj) 



where we define the following integral transform 



(3.28) 



E 



(/-«)! \il 



yM-1,0 



(-1)'-' /I! 
{l\f'-HLi{t) . 



M+l 



dS S* Gq 1 Q Q 



i+l(-uM-l 



(3.29) 



In the second step we have used the identity (|A.7p . which leads us back to the standard Laguerre 
polynomials. Taking into account the normalisation (|3.20p we arrive at the following final result, 



N-l 



1=0 



(3.30) 



which is proportional to the kernel of ordinary Laguerre polynomials p.l7p . The remaining kernel 
Hii can be expressed in terms of the two integral transforms which we have already computed, 



Hnit,s) = E T^i^rmris) - G--i\ L- ..,0 f ) 

1=0 K V V 



(3.31) 



This completes the computation of all (/c, /)-point correlation functions in the 2mm, together with eq. 

Although we will postpone the detailed analysis of the large- limit to future work let us mention 
the following nontrivial identity with respect to the kernel Hw. 



M-l.O 

0, A/-1 \ 0,...,0 



/=0 



(3.32) 



implying that limAr_5.oo Hn (t, s) = 0. Assuming that the sum converges and can be integrated piecewise 



this can be shown as follows. In appendix iBlwe verify that ^ (s) and x[*^'* {s) form a set of orthogonal 
functions with respect to the flat measure, see eq. p.34p . After multiplying both sides of eq. p.32p 
with p^l^'^^ (s) and integrating s over M_|_ we obtain 

which is consistent with eq. (]3.29p . 



dspf'\s)G, 



0,M~l \ 0,...,0 



(3.33) 
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The identity (I3.32P most likely implies that in the naive large- limit, meaning s, t and M fixed, 
the correlation function R^^^P factorises into s- and t-dependent parts regardless if Hqq vanishes or 
not since the determinant ()3.4p factorises into a k x k determinant incorporating Hqi and a I x I 
determinant comprising the block Hiq. Therefore the correlation functions of the singular values tj of 
the matrix Xi decouple and become the ones of the standard Wishart-Laguerre type. There may be 
other ways to obtain a nontrivial coupling between singular values Sa and tj in a more sophisticated 
large- limit (like in the so-called weak limit in |25|). 

We finally make contact again to the one- matrix model formulation (j2.12p . The following ortho- 
gonality relation which follows from eq. (j3.2p is explicitly verified in appendix [B] 



(3.34) 



in other words pri ^'^ (s) and (■s) constitute a set of biorthogonal functions with respect to a single 
variable with flat measure. In particular this relation results into the following property of the kernel 
Hqi, see eq. (|3.22p . that contains the two: 



/ ds'Hoi (s, s')Hoi (s', s") = Hoi (s, s") . 
Jo 

We have thus closed the circle back to the jpdf ()2.13p where we could directly replace 



(3.35) 



Am(s) det 

l<c.d<N 



A/,0 

0,M \ 0,...,0,d-l 



l<a,o<N l<c,d<N 



mi 



N-l 



-^-^ l<a,b<N 



(3.36) 



This uses the invariance property of determinants under addition of columns, and then proceeds 
with standard techniques to deduce the correlation functions. This directly leads from eq. (j2.13p to 
eq. (j3.10p for the jpdf. With the property (j3.35p of the kernel we deduce eq. (j3.26p from Dyson's 
theorem 



3.3 Spectral density, its moments and large-A^ scaling 

In this subsection we discuss in more detail the implications of our results for the spectral density of 
singular values. Starting from the expression eq. (I3.27P which we repeat here in two equivalent forms. 



Rf'\s) 



N-l I 

EE 

1=0 i,j=0 



(_l)i+i(/!)2 



M,0 



(Z - j)!(/ -i)!(i!)2(j!)Am ^o,m 
iFM(-j;l,...,l;s)Gf 



-lA'L 1 



j=0 



J'- 



, M+1 \ 0,...,0 



we can explicitly compute expectation values for the moments for finite-A^. Starting from the first 
expression we obtain 



-iy+'in)\i + j + ky.iu + k)\)''-^ 



N 



EE- 

Z=0 i,j=0 



(/-j)!(/-i)!(i!)2(j!)Am 



(3.37) 
(3.38) 
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Here we have normalised by eq. (I2.19P and we have used again the identity eq. (lA.Th for moments 
of the Meijer G-function. On the other hand using the compact expression for X^j^^\s) in the second 
formulation of the density we can obtain a more concise result in the following way: 



iA/, 1 



-y 



i=o ^ ^ j=o ^ -^^ 



1 ^V(fc+o!^^^+^ 



iV ^ V ?! 7 L 2nkV ' 1 - e-*"^ 



iV ^ A;! V l\ / V iV-/-l 

i=0 



(3.39) 



We use the convention that inverse powers of factorials of negative integers give zero, rather than 
using the Gamma- function everywhere. In the first step we have employed eq. (lA.Sp . Interestingly, 
the remaining sum can be further expressed in terms of a hyper geometric function if > A^, 

E[s'=] = (_i)A^-i^^i__L__i_ u-,2FM^x(k + 1, . . . , + 1, 1 - AT; 1, . . . , 1, _ AT + 1; 1) , (3.40) 

by extending the sum to infinity and comparing their Taylor series. Indeed this relation can be 
generalised to A; < A^. Notice that in this case the singular contributions in the hypergeometric 
function cancel with those in the Gamma- function in the denominator. 

In the ensuing discussion we will need in particular the first moment F^^^ for k = I when rescaling 



the density, which can be readily read off 

E[s] = Fjf^^ = A^^ . (3.41) 

It agrees with the known case for M = 1. An alternative short derivation for the first moment is 
sketched in appendix [Bj Higher moments easily follow from eq. (|3.39|) . e.g. for the second moment we 
have 

E[s2] = ^N^^ ((A + 1)^^+1 - (A - 1)^^+1) . (3.42) 

We illustrate our results for the density (13.27P by plotting it for various values of A and M. As 
a first example in fig. [T]the density is shown for M = 1, 2, 3 at fixed A = 4. Clearly it is mandatory 
to know the right scale dependence of the correlation functions on A and M. This means that after 
properly rescaling the bulk of the singular values is of order one, in order to be able to compare 
the density for finite A at different values of A^. In particular it is important to check the finite 
A-results against the limiting large-A behaviour for different M, which has been derived for products 
of quadratic [H [7] and rectangular matrices [H [9] . 

Let us explain our procedure. First we normalise our density to unity, using eq. (|2.19p . Then we 
rescale the density by its first moment, 

£,{M), ^ 1 (M) (M) f j^iM) ^ 



Rr>ix)^-F^-^Rr'[Fr>x) , (3.43) 
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Figure 1: Comparison of the density eq. (j3.27p Rl (s) for fixed = 4 without rescahng. The values 
M = 1,2,3 correspond to the top (blue), middle (red) and bottom (black) curve, respectively. 



so that the new density R^'' (x) has norm and first moment equal to unitjU. Notice that this rescaling 
is an alternative to the unfolding procedure onto the scale of the local mean level spacing. Instead 
of fixing the mean distance between two successive singular values, we fix here the singular values 
themselves, such that they are always of order one. This is exactly the macroscopic limit. 
Inserting eq. (|3.4ip we thus obtain a limiting density, which we denote by, 

p^^'\x)= lim i?f')(x)= lim N''-^r[^'Un^''x) , (3.44) 

that also has norm and first moment unity. Note that the known limiting density from the literature 
may still have to be rescaled accordingly. We can then compare p^^^\x) and r!i ^\x) for different 
values of at fixed M . 

This procedure can be illustrated with the simplest case M = 1, the well known Wishart Laguerre 
ensemble. At large-A^ we have 



-,(M=1). . 1 IAN -s 



hmi?r-^^(^)-^V^ ,G(4iV-.) (3.45) 

which is the A^-dependent Marchenko-Pastur density, and 0(x) denotes the Heaviside function. With 
the first moment being given by F^~^^ = N we thus have 



.(M=l). ^ 1- ^rO„(Af=l) . 1 /4 - X 



(x) = lim N^RY"-'' (Nx) = —\ 6(4 - x) , (3.46) 

7V->oo 27r V X 

for the rescaled and normalised density. This is the Marchenko-Pastur density with compact support 



IS 



on (0,4]. A comparison between p^^^^'^\x) and R\ ~ (x) for various values of = 3,4,5 and 10 
given in fig. [2] (left figure) . 

For M > 1 the limiting expression for the density is not as explicit as in eq. (I3.46p . Here we will 
follow the notation of [9] where a polynomial equation for the resolvent G{z) was derived, which we 
display for the case of quadratic matrices only: 

{z)) = z (z) - l) . (3.47) 



■^AU densities satisfying this property will be denotes with a hat 
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Figure 2: Left plot: the rescaled density '{Nx) eq. (f3:27|) for M = 1 and = 3,4,5,10 

(corresponding to red, black, cyan, and magenta, respectively) versus the limiting Marchenko-Pastur 
density eq. (^Mh (thick line). Right plot: the density eq. (lOTD for M = 2 rescaled as NR[^^"'^\n^s), 
for N = 3,4,5, 10 compared to the limiting density p^^^^'^\x) (thick line). 



The resolvent is related to the limiting spectral density by 



f^dxP^^ , (3.48) 
z — A 



Jo 

where z is outside the support of the spectral density. This relation can be inverted as follows: 

pW(A) = -i lim 9m f G^^^\X + ie)] . (3.49) 
vr e->0+ V / 

For M = 1 eq. (|3.47|) reduces to a quadratic equation, which after taking the discontinuity along the 
support eq. (|3.49p leads to eq. ()3.46p . without further rescaling necessary. 

Increasing to M = 2 the equation becomes cubic and we can still write out its solution, which is 
chosen subject to boundary conditions and to yield a real density on the support (0, 3^/2^]g 



G^''-'Hz) = -^{AZ'^\z)+A'J\z))=^{i-A^iz))^f'^ + A]/\z)) , 




'?,z 

A,(.) . I^yri-li^^j . (3.50) 

The density p(*^^^)(s) that is obtained from eq. (|3.50p by taking the discontinuity according to eq. 
(|3.49p (which happens to have the first moment equal to unity without further rescaling) is shown in 
fig. [2] (right figure) in comparison to our rescaled finite- result r\^^~'^\s) for various values of A^. 
As in the known case M = 1 we obtain a nice agreement for M = 2. An alternative derivation of 
^(M=2)^g^ via multiple orthogonal polynomials was quite recently presented in Ref. [31], with which 
our density agrees. 

As a final remark for larger M it might be useful to resolve the singularity of the density at the 
origin [9], liuis^o P^^Hs) ~ 5-^/(^+1)^ by changing variables, just as it is well known that for M = 1 



a change to squared variables maps the Marchenko-Pastur density eq. (I3.46P to the semi-circle. A 



*We thank Z. Burda for indicating how to determine the hmiting support. 
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Figure 3: The ergodic mutual information of multi-layered scattering MIMO channels with a fixed 
number of 2 clusters (M = 3) with a different number of scatters = 2, 4, 8. 

more detailed comparison to existing large- results for the density, its moments and support would 
require a careful asymptotic analysis of the special functions constituting our finite- density, and is 
postponed to future work. 

4 Application to telecommunication 

Consider a MIMO network with a single source and destination, both equipped with N antennas. 
Information transmitted by the source is conveyed to the destination via M — 1 successive clusters of 
scatterers, where each cluster (layer) is assumed to have N scattering objects. Such a channel model 
proposed in [8] is typical in modelling the indoor propagation of information between different floors 
[32]. 

We assume that the vector-valued transmitted signal propagates from the transmitter array to the 
first cluster, from the first to the second cluster, and so on, until it is received from the (M — l)-st 
cluster by the receiver antenna array. Each communication channel is described by a random complex 
Gaussian matrix, and as a result the effective channel of this multi-layered model equals the product 
matrix Pm, see eq. p.ip . 

For the described communication channels, the mutual information measured in units of the natural 
logarithm (nats) per second per Hertz is defined as 

N 

X(7) ^ Indet (Iat + j^PmPL) =J2^n[l + j^s^) , (4.1) 

a=l 

where 7 defines the average received signal-to- noise ratio per antenna which is a constant. We employ 
the distribution (|3.27p of squared singular values to compute its average. The quantity of interest 
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is called the ergodic mutual information of such channels. It is given by the expectation value of the 
random variable X(7) . Using the analogue of the expression (|3.37p from the previous section we have 



E[X(7)] = iVE|ln(l + ^s 



(_l)»+i(/!)2 
(_l)»+i(/!)2 



Jo 



ds Gn. AT . - 



s In 1 + 



EE 

1=0 ij=0 
N-l I 

^^^^J^^lW^^W^iW^^^'^^^ V o,o,j + i,...,j + i,i + j + i 



7 



■G 



M+2, 1 



0,1 



ATA/ 

7 



.(4.2) 



The last step is obtained by first replacing the logarithm by a Meijer G- function, eq. (IA.5p . and then 
applying the integral identity ()A.14p from the appendix. Note that the corresponding ergodic mutual 
information for the traditional MIMO channel model, i.e. M = 1, was derived in |33j . 

In order to get an independent confirmation of our analytical result (j4.2p we compare it to numerical 
simulations as follows. We plot the ergodic mutual information (j4.2|) in fig. [3] against Monte-Carlo 
simulations as a function of 7 in decibel (dB) for a given number of clusters, M — 1, as an example, 
with different numbers of scatters per cluster, N. Each simulated curve is obtained by averaging 
over 10^ independent realisations of Pm- The statistical error bar is smaller than the symbol for the 
simulation in our plots. The comparison in fig. [3] shows a 2-layered scattering channel, i.e. M = 3, 
with the number of scatters per layer varying from N = 2, 4, to 8. We have also compared our results 
to simulations for other values of A'^ and M. 



5 Conclusions and open questions 

In this paper we have derived the joint probability distribution function (jpdf) of singular values 
for any finite product of M quadratic random matrices of finite size N x N, with complex elements 
distributed according to a Gaussian distribution. This generalises the Wishart-Laguerre (also called 
chiral Gaussian) Unitary Ensemble which we recover for M = 1. Starting from the jpdf we have 
computed all fc-point density correlation functions of the singular values, by taking a detour over a 
two-matrix model like representation of the same model. In that way we showed that the jpdf being 
proportional to a Vandermonde times the determinant of Meijer G- functions represents a determinantal 
point process. Its kernel of orthogonal functions generates all /c-point functions in the standard way 
using Dyson's theorem. We also solved the auxiliary two-matrix model that couples a single matrix 
to the product of M matrices, by constructing the biorthogonal polynomials explicitly, as well as 
the corresponding four kernels with their integral transforms. On the way we found some nontrivial 
identities and integral representations of the Meijer G-function. 

The density of the singular values are discussed in more detail at finite N and M, including all 
its moments. We identified the macroscopic scaling to match the density with the known large-A^ 
results for the macroscopic density of singular values. As a further application we have computed the 
averaged mutual information for multi-layered scattering of MIMO channels and have compared them 
to Monte-Carlo simulations for small M and N. 

Previous results for the macroscopic large- density of the singular values of quadratic or rect- 
angular matrices and its expectation values of traces have mainly been obtained from probabilistic 
methods, in particular using free random variables. The explicit results that we have obtained for 
the jpdf and all correlation functions thereof open up the possibility of another direction. One can 
now investigate the microscopic scaling limits zooming into various parts of the spectrum, by per- 
forming the asymptotic analysis of the orthogonal polynomials and their integral transforms that we 
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computed. Since the ensemble represents a determinantal point process one can also investigate the 
limiting distribution of individual eigenvalues, e.g. by considering their Fredholm determinant repre- 
sentation. Moreover, one can now study the distributions of linear statistics such as the trace of the 
derived ensemble as well. 

Based on known results for the universality of the spectrum of random matrices we expect the 
following outcome of such an analysis. The bulk and the soft edge behaviour of the spectrum should be 
governed by the universal Sine- and Airy-kernel, respectively, after unfolding and scaling appropriately. 
This includes that at the soft edge we expect to find the Tracy- Widom distribution, and it will be very 
interesting to identify the right scaling for that. In contrast at the hard edge we expect to find new 
universality classes labelled by M. Already the way the macroscopic density diverges at the origin 
depends on it. This would fit into what was found recently for the complex eigenvalue spectrum of 
products of independent Ginibre matrices. 

Other open problems include a generalisation of our construction to products of rectangular matri- 
ces, which seems quite feasible. Also the inclusion of determinants (or characteristic polynomials) into 
the weight, which should be related to rectangular matrices, seems within reach. On the other hand 
going to non-Gaussian weight functions or investigating other symmetry classes with real or quater- 
nion real matrix elements are very challenging problems. The reason is that our method depends 
crucially on the Harish-Chandra-Itzykson-Zuber group integral to eliminate the angular variables. 
For the other symmetry classes no such explicit tool is presently at hand. 
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A Some integral identities for Meijer G-functions 



In this appendix we collect a few integral representations and identities for the so-called Meijer G- 
function. It is defined as 1231 



ai, 02, 
61, 62, ••• 



1 

27ri 



du 



m=™+ir(i-6,+n)m=„+ir(a, 



u 



u 



(A.l) 



The contour of integration C goes from — ioo to +ioo such that all poles of the Gamma functions 
related to the bj lie to the right of the path, and all poles related to the aj to the left of the path. We 
are in particular interested in the case 



m, 
0, m \ b 



X „ m 

n .7 — 1 



(A.2) 



Note that the function is symmetric in all its indices 61, ... ,6^- Special cases for small m are given 
by [23] 



•-^0,1 



2,0 

0,2 I 61,62 



2^(^.+62)/2^^^_^^(2Vi) , 



(A.3) 
(A.4) 
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and 



G 



1,2 
2,2 



1 1 

1 



ln(l + z) , 



(A.5) 



where K is the modified Bessel function of the second kind equivalently known as the Macdonald 
function. It is possible to absorb powers of the argument of the Meijer G-function, due to the following 
shift 



(^m,n I 0-1 + ^, • • • , Qp -\- k 

P'" \ h + k, b„ + k 



(A.6) 



The integral identities we collect here are for the moments of the Meijer G-function: 



j^^n-l^m, 



r n^M, 1 / -J 
Ub b \ 0,...,0 



n 



n + 

i=i 

'n])^r{j -n + e) 
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£^0 



r(-n + e) 



i-iy 



(n!) 



r(n - j + 1) 



(A.7) 



(A. 



We only state two particular cases, for the most general setting see e.g. [34]. In the last step of the 
second identity we have employed the recursion of the Gamma-function, T{j — n-\- e) = 

r(-n + e)nCo(^ + 

The second integral identity needed in section [2] follows easily from the representation eq. (1A.2P : 



V \v 



-s/v^m, 
^ ^0,m \ bi,...,b„ 



c ^-^ 



G; 



m+l, 



— [ duf\ T{bi - u) s"r((i -1-u) 
J 



c J=i 



0, m+l I bi,...,bm.,d~l 



(A.9) 



In the first line we simply substituted t ^ v = s/t to obtain a second version of the identity. The 
rest follows from the definition of the Meijer G- and the Gamma-function. Notice that for d = 1 and 
bi = . . . = bm = this recursion for the Meijer G-function was already derived in jl6j . 

We are now prepared to show the multiple integral representation of the Meijer G-function that 
we need in the derivation of the jpdf in section [21 It slightly generalises the representation found in 
1161. The statement is that 



g: 



m, 

0,m I 0,...,0,f) 



Xo 



^ dX2 



dxi ( Xi \ I ax2 

Xi \XoJ Jq X2 



f^m—l 



Xm-1 



j/Xj-l 



(A.IO) 



for m > 1. Comparing this equation to eq. (|2.11|) with b = d — 1, we work with m = M squared 
singular values, Xj = (Ac"'^)^, where we introduce a dummy variable xq. The variable xq will be set to 
unity for our original purposes. However, it will be useful when applying the identity to eq. (I2.2ip . 
Our proof goes by induction in m. For m = 2 [23] we have 



dxi ( xiY _fi _f2 
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(A.ll) 
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where the last step is due to eq. (IA.4h . This leads to 

-■^ dx 



L 
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2^2,0 
•-^0,2 I 0,b 



xz/x2 /^3, 

~ '-'0,3 I 0,6,0 



X3 

Xq 



(A.12) 



for m = 3, where we have apphed the identity eq. ()A.9P in its second form shown in the first Une, 
with d = 1 and v = x^jx^. For the induction step m — 1 — t- m we simply have to repeat the same 
procedure, which follows easily from the very same identity 



^ 



m i^m, 



0,m I 0,...,0,6 



m+1,0 

0,m+l I 0,...,0,6 



^m+1 



(A.13) 



which completes the proof. 

Note that the same identity (jA.lOp can be used to provide the second step in eq. (j2.2ip . when 
setting 6 = 0, m = M— 1 and shifting the indices of the variables = (A^''^)^ for j = 1, . . . ,M. 
This is the reason why xq is useful. 

Finally we state an integral identity needed in section U concerning the integral of two Meijer 
G- functions, 



ds G 



1,2 
2,2 



1 1 
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7 \ pM,0 



J, 



-1,0 



N 



(A.14) 

Notice that it is a particular choice of a general formula |34] . In order to arrive at eq. (|4.2|) we apply 
the shift (IXell . 



B Orthogonality check and first moment 

In this appendix we explicitly confirm both the orthogonality ()3.2p of the bOP Pri ^\s) and ql^^\t) 

with respect to two variables, as well as the fact that pi^^\s) and x[^^\-^) constitute a set of biorthog- 
onal functions with respect to one variable with fiat measure, eq. p.34p . Although being true by 
construction we will see the orthogonality ultimately boils down to the standard orthogonality of 
Laguerre polynomials. 

The biorthogonal polynomials that were constructed in subsection 13.11 using the bimoment matrix 
must automatically satisfy the orthogonality relation p.2p . We will check this here independently, 
which implies at the same time that one of the polynomials and the integral transform (j3.23p of the 
other are orthogonal functions (as they should be, in order to constitute proper kernels): 

POO POO 

J ds pf^ is) xf ) (.) = j dsdt W t) pf^ {s)qf'^ it) 

= {^l)''+'S.,. (B.l) 

Here we have used the identity ()A.7p for moments of the Meijer G-function, which cancels the extra 
factorials in the generalised Laguerre polynomials p\^\s), see eq. (j3.2ip after the first integration 
over s' = s/t. The last step follows from the known orthogonality of Laguerre type. 
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In the second part of this appendix we provide a much simpler, probabiHstic argument that leads 
to the first moment eq. (j3.4ip . Noting that 



E 



NIm , for j = l,...,L 



(B.2) 



we have that 



1e 

N 



N 



N 



a=l 
1 / 



1 



-iV^"Tr(/^)=iV^ , 



(B.3) 



where we have reordered successively the Xj under the trace in the first line. 
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